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SUMMARY 

Two finite element analysis codes are used to model the effects 
of cooling on the dimensional stability of formed and bonded 
composite parts. The two analysis routines, one h-version and one 
p-version, are compared for modeling time, analysis execution time, 
and exactness of solution as compared to actual test results. A 
recommended procedure for predicting temperature effects on 
composite parts is presented, based on the results of this study. 

INTRODUCTION 

McDonnell Aircraft Company is actively involved in the research 
and development of advanced composite structures. The goals of 
lighter weight, lower cost, and more survivable composite aircraft 
structures are being pursued by several research projects. 

Innovative designs are being studied in conjunction with emerging 
materials to realize these goals. 

Several studies utilizing the thermoplastic composite material 
system are currently in progress. Design, manufacturing, and 
testing of sub-scale and full-scale components are part of the 
research activities. An interesting phenomenon which has arisen 
during fabrication and bonding is a dimensional warpage or 
“sprlngback" of the sub and full scale thermoplastic components. 

This warpage is caused by the difference in the coefficients of 
thermal expansion (CTE) in the principal directions of the 
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the rmop 1 as t i c material. As the part is cooled, the difference in 
the in-plane and through thickness CTE's causes a build-up of 
thermal stresses. The residual thermal stresses cause a warpage 
when the part is removed from the tool. In some instances, the part 
deflected off the tool at the edges as it was cooled. The resultant 
sprlngback can occur in thermoplastic, thermoset, or any other 
composite material with different directional CTE's. 

Related investigations have been completed to attempt to 
predict the amount of springback in formed or bonded composite 
parts. Material properties have been shown to be important in the 
prediction of springback magnitudes, references [1] and [2]. The 
through thickness coefficient of thermal expansion is an influential 
material property. Fahmy and Ragai-Ellozy derived a formula to 
estimate the through thickness CTE in reference [3]. The formula i 
based on known in-plane thermal properties. 

Finite element solutions have h;en used to predict springback 
for simple geometries with some appirent success, references 
[4]-[6]. The majority of work has demonstrated that simple 
geometries can be modeled for springback with a reasonable degree of 
accuracy. However, more complex shapes are much more difficult to 
model . 


Additional work is being completed in research activities at 
McDonnell Aircraft Company. The degree of springback has been shown 
to be linked to the build-up of residual thermal stresses. The 
effects of tooling temperature and ply orientation on the magnitude 
of these residual stresses are being investigated. The residual 
stresses have been shown to dramatically affect the interlaminar 
tensile strength of thermoplastic test specimens. 

PROBLEM STATEMENT 

The need exists to develop a reliable technique to predict 
springback in formed and bonded composite structures. Current 
procedures for springback compensation involve the reworking of 
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tools after the magnitude of springback is observed on a final part. 
This operation is costly and must often be repeated several times, 
depending on the complexity of the part. Tooling can be designed to 
compensate for springback if a reliable analysis method can be 
formulated. The compensated tool would then be capable of producing 
parts which are dimensionally correct the first time and eliminate 
the need for rework. 

The use of finite element methods to predict springback has 
shown promise. H-version finite element codes have been used to 
model simple geometries with reasonable agreement with actual data. 
Lesser degrees of accuracy and reliability are realized for more 
complex geometries. The p-version finite element approach is an 
alternative to the current h-version methods. The p-version 
formulation increases the degree of polynomial shape function for 
each element as the mesh is refined. The number of elements remains 
constant and, therefore, requires less modeling time than an 
h-version analysis. Thus far, however, very little work has been 
done using p-version finite element codes to predict springback 
magnitudes. 

This paper compares the h and p-version finite element codes 
for predicting springback magnitudes. Simple and complex geometries 
of Isotropic and orthotropic material properties are modeled. The 
two analytical methods are then compared for modeling time, analysis 
execution time, and exactness of solution as compared to actual test 
resul ts. 


DESCRIPTION OF ANALYSIS CODES 

The two finite element codes chosen for comparison are COSTAR 
(h-version) and MSC/PROBE (p-version). 

COSTAR, the Composite STructural Analysis Routine, was 
developed by Jon Goering of the Structural Research Department of 
McDonnell Aircraft Company and is widely used in Research and 
Development projects. COSTAR elements are assembled and formulated 
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In a method similar to MSC/NASTRAN. Beam, shell, plate, and solid 
elements are available in linear, quadratic, and cubic forms. 

COSTAR Is particularly well suited for the analysis of composite 
structures. The program has the ability to consider laminated and 
fully three-dimensional reinforced composites, as well as Isotropic 
and orthotropic materials, reference [7]. The majority of figures 
In this paper were generated by the COSTAR post-processor. 

MSC / PROBE uses the p-version approach to the finite element 
solution by increasing the degree of polynomial shape function of 
the element to the point of solution convergence. This feature 
often decreases modeling and solution time due to the requirement of 
fewer elements. 

MSC/ PROBE offers an error estimation routine in which the 
finite element solution is compared to the exact problem solution. 
The exact solution Is based on the Principle of Virtual Work and is 
of the form; 


( 1 ) 

In which is a bilinear expression for the virtual work of 

Internal stresses and is a linear functional of the virtua’ 

work of external forces. Error estimation is based on solving t. 
equation such that the difference between the finite element 
solution and the exact solution, based on the Principle of Virtual 
Work, is minimized. This can be expressed as; 


( 2 ) 

in which is the exact solution of Equation (1) and is the 

user defined finite element solution. An excellent source of 
Information on the formulations and methodology used in MSC/PROBE 
can be found in reference [8], 

Other features such as the ability to calculate stress 
Intensity factors around singularity points are available. The 
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author would refer readers to the MSC/PROBE User's Manual for 
detailed descriptions of other capabilities, reference [9]. 


BASELINE COMPARISONS 

CASE I: SIMPLE GEOMETRY - ISOTROPIC MATERIAL 

The Initial test of the ability of COSTAR and MSC/PROBE to 
model sprlngback was a simple isotropic hat section. Figure 1. The 
material properties used are as follows: 


E = 0.5 x 10 6 psi NU = 0.35 

CTE = 0.31 x 10 -6 1 / 0 F 


The model discretization required for the two programs can be seen 
in Figures 2a and 2b. The COSTAR mo-'al required 388 elements and 
980 degrees of freedom while the M'3/PR0BE model required 32 
elements and 1048 degrees of freedom at a p-level of 5. The applied 
load was a -500°F temperature gradient to represent the cooldown of 
the part from consolidation temperature. The constrained model is 
shown In Figure 1. These load and constraint conditions were used 
on all models for the duration of the study. 

The displacement results for the simple geometry, isotropic 
case can be seen in Table 1. As can be seen, COSTAR and MSC/PROBE 
give the same displacement results for this configuration. The 
error estimation for MSC/PROBE p-levels 1-8 is seen in Table 2. The 
execution time for the two analyses is very similar for the single 
COSTAR run and the MSC/PROBE run for p = 5 only. The final 
deflected geometry can be seen in Figure 3. 

CASE II: SIMPLE GEOMETRY - ORTHOTROPIC MATERIAL 

The effect of material property on the results of the two 
analysis programs was then studied. The isotropic material of Case 
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I was changed to the following orthotropic properties: 


Eji = 2.24 x 10 7 psi 
e 22 “ 1*30 x 10® psi 
E 33 - 1.30 x 10® psi 


NU 12 = 0.35 
NU 2 3 = 0.48 
NU 31 = 0.02 


CTE n = 0.10 x 10"® 1 / °F G = 0.79 x 10 6 psi 

CTE 2 2 - 0.165 x 10" 4 1 / 0 F 
CTE 33 - 0.175 x 10" 4 1 / 0 F 


All other model parameters were identical to Case I. The results 
again show very close agreement between the two models. Table 3. 
Error estimations for the Case II r.odel are s.milar to the Casr J 
calculations seen in Table 2. The deflected geometry can be seen 

Figure 4. 


CASE III: REFINED MESH - ORTHOTROPIC MATERIAL 

The effect of mesh discretization on the results of the 
MSC/PROBE model was demonstrated by performing an h-p extension at 
the two radii. The number of elements was increased to 40 which can 
be seen In Figure 5. The results of this analysis again show dost 
agreement between the two programs. Table 4 shows the displacements 
calculated from this analysis which are nearly identical to the 
results of Case II, Table 3. The effect of the h-p extension on the 
energy norm error can be seen in Table 5. The calculated error 
between the exact solution and finite element solution is decreased 
for the same p-levels of Table 2. 

CASES I - III RESULTS SUMMARY 

The results of Cases I - III show nearly exact agreement 
between COSTAR and MSC/PROBE in the calculation of temperature 
Induced displacements. Analysis execution times between a single 

COSTAR run and an MSC/PROBE p-level 5 run are similar. MSC/PROBE 

1544 



provides an error estimation function which can be used to evaluate 
the accuracy of the MSC/PROBE model. An h-p extension on the 
MSC/ PROBE model was shown to decrease the percent error in the 
energy norm for the simple hat geometry. 

STIFFENED SKIN ANALYSIS 


The results of the previous case studies indicate that COSTAR 
and MSC/PROBE will predict similar springback magnitudes for a 
simple geometry. The current research problem, however, is the 
springback encountered in the bonding of a thermoplastic 
hat-stiffened skin. This more complex geometry is modeled in an 
attempt to accurately predict the expected amount of springback. 

Material properties for the corrugated panel and the flat skin 
panel to be bonded were determined using the program THICKLAM, 
reference 10. THICKLAM utilizes Classical Lamination Theory to 
derive three-dimensional elastic constants of composite laminates 
using the formulae of Jones and Sun, references 11 and 12. The 
three-dimensional properties are necessary due to the cross 
sectional geometry of the model. The model X axis follows a 
material 90° ply and the model Y axis is through the thickness of 
the skin. Therefore El of the model is Ey of the material and E2 of 
the model is Ez of the material, Figure 6. The material properties 
used for the corrugation and flat skin panels are shown in Tables 6a 
and 6b. 


The COSTAR and MSC/PROBE finite element models are shown in 
Figures 7a and 7b. The COSTAR model required approximately 1500 
elements and 2325 degrees of freedom while the MSC/PROBE model 
required 153 elements and 3400 degrees of freedom for p = 4. The 
material properties of Table 6 were used for the corresponding parts 
of the model and the proper local material properties were input 
Into the MSC/PROBE model. The applied load was a -500°F temperature 
gradient as In the previous case studies. The constraints are again 
similar to the previous studies and are shown in Figure 6. 
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The results of this study again show close agreement between 
the two finite element programs, even for the more complex geometry. 
Table 7 shows the deflections for the two models. Figure 8 shows 
the deflected geometry of the COSTAR model which is essentially the 
same as the MSC/PROBE model. The more sophisticated model produced 
a greater error as can be seen in the calculations of Table 8. The 
COSTAR model required approximately one minute of execution time for 
the solution. The MSC/PROBE model required approximately forty 
minutes of execution time for all p-levels 1-8. 

SUMMARY OF RESULTS AND CONCLUSIONS 

The COSTAR and MSC/PROBE finite element programs predict almost 
exactly the same magnitudes of springback for a range of geometries 
and materials. Less time was required to construct the MSC/PROBE 
models due to the fewer number of elements. The addition of the 
orthotropic material angle in MSC/PROBE Version 4.1 further reduces 
the amount of modeling time required. 

A single execution of the COSTAR model was used throughout this 
project as a comparison to the eight p-levels solutions available 
with MSC/PROBE. This is an important factor due to the error 
encountered in the final stiffened skin model. The COSTAR model 
must be re-meshed in order to improve the solution. The MSC/PROBE 
model would require less extensive alterations and, therefore, 
requires less time to improve the solution. 

The two programs show good agreement with actual results. A 
stiffened panel of the same dimensions and material properties as 
those modeled was previously fabricated and bonded. The total 
springback measured for the panel was 0.29 inches. The magnitude of 
the springback predicted by both COSTAR and MSC/PROBE is 0.22 inches 
for a difference of 24%. 


FUTURE WORK 

The finite element methods used for this study have been shown 
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to predict the correct deflected geometry of a thermoplastic 
stiffened skin. Current research has shown that such factors as 
material properties and geometry have a great impact on the 
magnitude of springback in a composite part. Other physical and 
model parameters should be investigated in order to improve the 
predicted magnitudes of springback. 

The bondline of the stiffened skin was thought to have little 
effect on springback and was not modeled. This may be an important 
omission from the model. Material properties, particularly the 
through thickness CTE, must be accurately determined for use in 
springback prediction. Thermoset materials experience a chemical 
shrinkage during cure which effectively increases the through 
thickness CTE. This effect must be accurately quantified if 
thermoset springback prediction is to be accomplished. 

Geometry plays an important role in springback prediction. 
Current models require two dimensional representation only. 
Eventually this must be extended to three dimensions. Plane strain 
plate elements are used for the two dimensional cases and seem to be 
the most reliable. Other element formulations may prove to be more 
accurate. The type of element to be used for the three dimensional 
case has not been investigated. 

The solution used for this project was for an applied load of a 
single temperature step. Material properties were held constant 
over the entire interval. The actual case would change the material 
properties as the part is cooled from the stress free state. A 
piecewise linear solution may be more realistic than the single 
Interval used for this study. Several temperature steps should be 
applied to the model with material properties varied accordingly for 
each interval. This load and material property condition would more 
realistically model the actual behavior of the part as it cools. 

RECOMMENDED APPROACH FOR SPRINGBACK PREDICTIONS 


1. Verify all geometry and use on the finite element model as 
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accurately as possible. 


2. Components with three-dimensional curvature should be sectioned 
Into appropriate two-dimensional sections, if at all possible. 

Solid element modeling of complex three-dimensional parts is 
possible, but the advantages over 2-D modeling are limited and 
normally not worth the additional modeling time. 

3. Two-dimensional generalized plane strain elements are 
recommended. These elements seem to most closely model the actual 
stress/strain behavior of the material, particularly if the material 
Is a composite. 

4. At least five elements should be used to define a radius. Four 
elements through the thickness are necessary for an h-version finite 
element code while only one element is required for a p-version 
code. 

5. Temperature dependent orthotropic properties should be used to 
define the materials. The properties should be broken into 
corresponding temperature intervals, the size of which will depend 
on the material and overall change in temperature. 

6. The finite element model is then run for each of the defined 
temperature intervals with corresponding mtcerlal properties. The 
final solution is the summation of the resulting deflections from 
each analysis run. 
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TABLE 1. DISPLACEMENT RESULTS FOR CASE I, ISOTROPIC MATERIAL 


MSC/PROBE 

COSTAR 

Node# 

Delta X 

Delta Y 

Node# 

Delta X 

Delta Y 

1 

0.0 

-0.0200 

1 

0.0 

-0.01999 

5 

0.0 

-0.212 

5 

0.0 

-0.02118 

931 

-0.0128 

-0.00166 

437 

-0.01288 

-0.00156 

936 

-0.137 

-0.00242 

141 

-0.01378 

-0.00234 

911 

-0.0502 

0.0 

911 

-0.05022 

0.0 

915 

-0.0502 

-0.00120 

915 

-0.05022 

-0.00197 


GP1 4-0169-1 1-D/cm 


TABLE 2. ERROR ESTIMATION FOR CASE 1, ISOTROPIC MATERIAL 


P 

Global 

DOF 

Delta Y 

Extrapolated 

Energy 

Convergence 

Rate 

Percent Rel 
Error Est'd 

1 

120 

0.102792Q-01 

0.000000D+00 

0.00 


2 

304 

0.4901 44D-05 

0.000000D+00 

9.65 

60.25 

3 

488 

0.431 920D-05 

0.431 51 5D-05 

1.25 

44.84 

4 ! 

736 

0.408346D-05 

0.385223D-05 

0.96 

36.81 

5 1 

1048 

0.401 905D-05 

0.3984160-05 

0.40 

34.29 

6 

1424 

0.381 051 D-05 

0.356409D-05 

2.22 

24.41 

7 

1864 

0.377486D— 05 

0.376470D-05 

0.68 

22.29 

8 

2368 

0.370690D-05 

0.35961 9D-05 

2.00 

17.55 


G PI 4-01 69-1 2-D/cm 
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TABLE 3. DISPLACEMENT RESULTS FOR CASE II, ORTHOTROPIC MATERIAL 


MSC/PROBE 

COSTAR 

Node# 

Delta X 

Delta Y 

Node# 

Delta X 

Delta Y 

1 

0.0 

0.00371 

1 

0.0 

0.00372 

5 

0.0 

0.00300 

5 

0.0 

-0.00301 

931 

-0.0143 

0.00082 

437 

-0.01443 

-0.00078 

936 

-0.0144 

-0.00010 

441 

-0.01457 

-0.00012 

911 

-0.01 52 

0.0 

911 

-0.01533 

0.0 

915 

-0.0152 

-0.00071 

915 

-0.01533 

-0.00071 


GP1 4-0 169-1 3-D/cm 


TABLE 4. DISPLACEMENT RESULTS FOR 
CASE III, ORTHOTROPIC MATERIAL, 
REVISED MSC/PROBE MESH 


MSC/PROBE 

Node # 

Delta X 

Delta Y 

1 

0.0 

0.00372 

5 

0.0 

0.00301 

931 

-0.0148 

0.00034 

936 

-0.0148 

-0.00045 

911 

-0.0153 

0.0 

915 

-0.0153 

-0.00071 


GP14-0169-14-D/cm 
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TABLE 5. ERROR ESTIMATION FOR CASE III, REVISED MSC/PROBE MODEL, 

ORTHOTROPIC MATERIAL 


p 

Global 

DOF 

Total 

Energy 

Extrapolated 

Energy 

Convergence 

Rate 

Percent Rel 
Error Est'd 

1 

152 

0.682723D-02 

O.OOOOOOD+OO 

0.00 


2 

348 

0.407940D-05 

O.OOOOOOD+OO 

9.82 

47.89 

3 

616 

0.391 884D-05 

0.391814D-05 

0.50 

42.54 

4 

928 

0.370865D-05 

0.354308D-05 

1.05 

34.29 

5 

1320 

0.355625D-05 

0.266208D-05 

1.41 

26.77 

6 

1792 

0.344546D-05 

0.282333D-05 

2.05 

19.57 

7 

2344 

0.338463D-05 

0.327377D-05 

2.43 

14.12 

8 

2976 

0.335552D-05 

0.331 843D-05 

2.43 

10.57 


GP14-0169-15-D/rr 


TABLE 6a. MATERIAL PROPERTY CALCULATIONS FOR THE CORRUGATION 


THICKLAM assumes that the composite is a balanced, symmetric laminate of one material system. 


Input 


Material System 

IM7/PEEK 

Lamina Orientation 


Angle 1 

0° 

Angle 2 

45° 

Angle 3 

-45° 

Angle 4 

90° 

Angle 5 


Angles 


Angle 7 


Angle S 


Lamina Data 


Ei 

2.24E+07 

e 2 

1 .30E+06 

e 3 

1.30E+06 

Pl2 

0.350 

Pl3 

0.020 

M-23 

0.480 

G i2 

7.90E+05 

Gl3 

7.90E+05 

G 23 

4.50E+05 

^21 

0.020 

^31 

0.001 

^32 

0.480 

A 

2.01308E-? 0 


Temperature 

RT 

Moisture Content 

Dry 

Percent Thickness 
Angle 1 
Angle 2 
Angle 3 
Angle 4 

30.9% 

30.8% 

30.8% 

7.7% 

Angle 5 
Angle 6 
Angle 7 
Jllfl Angle 8 



Total Thickness 

0.07 


Output 


Lamina (C) Matrix 


C11 

2.26E+07 

C12 

6.14E+05 

C13 

3.21 E+05 

C22 

1.71E+06 

C23 

8.20E+05 

C33 

1.69E+06 

C44 

4.50E+05 

C55 

7.90E+05 

C66 

7.90E+05 

Laminate Data 


Ex (E 3 ) 

9.46E+06 

Ey (E,) 

5.38E+06 

Ez (E 2 ) 

1.63E+06 

Pxy (P3l) 

0.541 

Pxz (jJ-32) 

0.102 

Pyz (Pi 2 ) 

0.278 

G xy (G 3 i) 

3.86E+06 

Gxz (G 32 ) 

6.12E+05 

Gyz (G 12 ) 

5.39E+05 

Pyx (Pl3) 

0.307 

Pzx (P23) 

0.018 

Pzy (P 2 l) 

0.084 


GP1 4-01 69-1 6-D/rr 
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TABLE 6b. MATERIAL PROPERTY CALCULATIONS FOR THE FLAT SKIN PANEL 

THICKLAM assumes that the composite is a balanced, symmetric laminate of one material system. 

Input Output 


Material System 

IM7/PEEK 

Lamina Orientation 


Angle 1 

0° 

Angle 2 

45° 

Angle 3 

-45° 

Angle 4 

90° 

Angle 


Angle 8 ii-HH 

Annb T •••••: 


/M iy It7 ; ; 

Angle 8 


Lamina Data 


Ei 

2.24E+07 

e 2 

1.30E+06 

e 3 

1.30E+06 

Pi 2 

0.350 

Pi 3 

0.020 

P23 

0.480 

Gi 2 

7.90E+05 

G 13 

7.90E+05 

^23 1 

4.50E+05 

P21 

0.020 

P31 

0.001 

P32 

0.480 

A 

2.01308E-20 


Temperature 

RT 

Moisture Content 

Dry 

Percent Thickness 
Angle 1 
Angle 2 
Angle 3 
Angle 4 

30.9% 

30.8% 

30.8% 

7.7% 

Angte 
i Angle® 
Angle? 

:||||:Ahgle^8:|pp| 

■ 

Total Thickness 

0.07 


Lamina (C) Matrix 


on 

2.26E+07 

Cl 2 

6.14E+05 

C13 

3.21 E+05 

C22 

1.71E+06 

C23 

8.20E+05 

C33 

1.69E+06 

C44 

4.50E+05 

C55 

7.90E+05 

C66 

7.90E+05 

Laminate Data 


Ex (E 3 ) 

1.14E+07 

E y (E,) 

7.37E+06 

Ez (E 2 ) 

1 .64E+06 

Pxy (P31) 

0.308 

Pxz (P32) 

0.194 

Pyz (P 12 ) 

0.305 

Gxy (G3l) 

2.78E+06 

Gxz (G 32 ) 

6.07E+05 

Gyz (G 12 ) 

5.44E+05 

Pyx (P13) 

0.200 

Pzx (H23) 

0.028 

Pzy (P 21 ) 

0.068 


TABLE 7. DISPLACEMENT RESULTS FOR THE FULL STIFFENED SKIN 


MSC/PROBE 

COSTAR 

Node# 

Delta X 

Delta Y 

Node# 

Delta X 

Delta Y 

236 

-0.00952 

0.0949 

1 

-0.00728 

0.10900 

229 

-0.01100 

0.0940 

11 

-0.01196 

0.10660 

363 

-0.00580 

0.0616 

103 

-0.00564 

0.06158 

237 

-0.00793 

0.0603 

116 

-0.00637 

0.05993 


GP1 4-01 69-1 8-D/an 


TABLE 8. ERROR ESTIMATION FOR THE HAT-STIFFENED SKIN 


P 

Global 

DOF 

Total 

Energy 

Extrapolated 

Energy 

Convergence 

Rate 

Percent Rel 
Error Est'd 

4 

3467 

0.171003D+00 

O.OOOOOOD+OO 

0.00 

14.28 

5 

4947 

0.169596D+00 

0. 000000 D+00 

1.49 

10.95 

6 

6733 

0.168854D+00 

0.167585D+00 

1.49 

8.70 


GP14-0169-19-D/rr 


1554 























• •• 



G P 1 4-0 1 69-4-0/kmg 


Figure 3. Deflected Geometry for Case ! ' 
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Figure 4. Deflected Geometry for Case II, Orthotropic Material 
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Figure 5. MSC/PROBE Model for Case III, Refined Mesh 
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Figure 8. Deflected Geometry for the Hat-Stiffened Skin 
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